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INTRODUCTION 


This  report  is  concerned  with  the  use  of  variational  principles  to  solve 
a  mixed  boundary  and  initial  value  problem.  From  a  previous  paper!  we  under¬ 
stand  that  the  far  end  conditions  are  not  imposed  for  solutions  in  an  initial 
value  problem.  This  implies  that  the  boundary  value  problem  can  be  solved  in 
strips  of  arbitrarily  chosen  intervals  of  time.  The  size  of  the  computation 
can  be  reduced  substantially  if  the  time  interval  taken  is  sufficiently  large 
and  number  of  strips  small.  This  depends,  of  course,  on  the  accuracy  of  the 
method . 

A  procedure  can  be  obtained  for  a  recursive  relationship  in  the  time 
domain,  where  the  final  conditions  of  first  strip  can  be  regarded  as  initial 
conditions  of  the  second  strip.  These  recursive  solutions  can  be  obtained  by 
using  variational  principles  with  the  aid  of  the  bicubic  Hermite  polynomial 
spline  functions  as  finite  elements. 

ESTIMATION 

A  dynamical  system  can  be  modeled  by  the  following  partial  differential 
equation. 

L(C)  ya(c)  =  -Q(G)  (1) 

with  appropriate  boundary  and  initial  conditions.  In  the  above  equation  L  is 
a  linear  operator,  in  both  spatial  and  temporal  domain,  ya  is  the  dependent 
variable,  Q  is  a  forcing  function  and  £  represents  all  Independent  variables, 

Ishen,  C.  N.  and  Wu,  J.  J.,  "A  New  Variational  Method  for  Initial  Value 
Problems  Using  Piecewise  Hermite  Polynomial  Spline  Functions,"  presented  at 
the  1981  Army  Numerical  Analysis  and  Computer  Conference,  Huntsville,  AL, 
February  1981. 
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both  spatial  and  temporal. 

The  Inner  product  <  >  of  an  adjoint  forcing  function  Q  and  the  solution 
(yaU))  of  Eq.  (1)  can  be  used  for  the  purpose  of  estimation.  This  inner 
product  is 

G[ya]  =  <Qya>  (2) 

The  estimate  y,  which  differs  from  the  actual  solution  ya  of  Eq.  (1)  by  an 
increment 

<Sy  =  y  -  ya  »  (3) 

then  becomes 

G[y]  =  G[ya  +  Sy] 

*  <Q>  (ya+(Sy)> 

=  G[ya]  +  <Q,<Sy>  (4) 

which  is  in  error  to  first  order  in  6y  and  Q.  This  is  undesirable  because  the 
error  depends  on  the  variation  6y  which  is  supposed  to  be  arbitrary.  Thus  the 
estimate  will  not  be  accurate. 

THE  VARIATIONAL  PRINCIPLE 

A  more  accurate  estimate  can  be  made  by  constructing  a  variational 
principle1  for  Eq.  (2).  By  using  the  adjoint  variable  y  as  a  Lagrange 
multiplier  for  Eq.  (1)  added  to  G[y]  we  have 

J[y,y]  =  G[y]  +  <y,(QfLy)> 

=  <Q,y>  +  <y,Q>  +  <y,Ly>  (5) 


1Shen,  C.  N.  and  Wu,  J.  J.,  "A  New  Variational  Method  for  Initial  Value 
Problems  Using  Piecewise  Hermite  Polynomial  Spline  Functions,”  presented  at 
the  1981  Army  Numerical  Analysis  and  Computer  Conference,  Huntsville,  AL, 
February  1981. 
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In  order  that  J  be  a  variational  principle  for  G  the  following  requirements 
must  be  satisfied. 

(a)  J  is  stationary  about  the  function  ys  which  satisfies  the  relation  in 
Eq.  (1). 

Lys  =  “Q  (6) 

(b)  The  stationary  value  of  J  deduced  from  Eqs.  (2)  through  (5)  is 

Jly.y]  =  G[ys]  ->■  G[ya]  (7) 

Consider  first  the  stationarity  of  J  by  taking  the  variation  of  Eq.  (5) 

<5J  =  <Q,  <5y>  +  <<5y,Q>  +  <<5y,Ly>  +  <y,L<5y> 

=  <6y,(Ly+Q)>  +  <Sy,(Ly+Q)> 

-  <6y,Ly>  +  <y,L6y>  (8) 

We  will  make  an  effort  later  to  impose  certain  conditions  in  order  that 
the  following  equality  holds: 

<y,L6y>  =  <6y,Ly>  (9) 

where  L  is  the  adjoint  operator. 

By  combining  Eqs.  (8)  and  (9)  one  obtains 

6 J  =  <6y,(Ly+Q)>  +  <6y,(Ly+Q)>  =  0  (10) 

Since  the  variations  6y  and  6y  are  arbitrary  it  leads  to  the  requirement  that 
the  stationary  values  ys  and  ys  must  satisfy 


O' 

1 

1! 

W 

(ll) 

Lys  =  -Q 

(12) 

Equation  (11)  is  the  same  as  Eq.  (6),  therefore  J  is  stationary  about  the 
function  ys. 
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Equation  (12)  is  the  adjoint  equation  in  terms  of  the  adjoint  operator  L, 
the  adjoint  variable  y,  and  the  adjoint  forcing  function  Q. 

It  is  noted  that  6J  in  Eq.  (10)  vanishes  and  is  independent  of  the 
arbitrary  variations  6y  and  Sy,  in  contrast  with  Eq.  (4) ,  where  SG  is  in  erroi 
to  the  first  order  in  Sy.  By  using  SJ  instead  of  SG  one  can  claim  that  the 
estimate  is  more  accurate  and  free  from  the  arbitrary  variations. 

Using  the  relationship  in  Eq.  (11)  the  stationary  value  of  J  from  Eq.  (5) 

is 

J[ys>ysl  =  <Q>ys>  +  <ys>Q>  +  <ys>Lys>  =  o[ys]  d3) 

Since  J  is  stationary  and  SJ  0,  then 

G[ys]  +  G[ya]  (14) 

which  is  the  requirement  given  in  Eq.  (7). 

It  is  noted  that  Eq.  (10)  contains  no  boundary  terms  to  be  satisfied. 

This  bears  an  important  point  in  the  future  discussion. 

BILINEAR  CONCOMITANT 

We  will  find  out  the  conditions  for  the  assumed  equality  in  Eq.  (9)  to  be 

o 

true.  Let  us  consider  the  following  bilinear  concomitant : * 

D  =  <y ,Ly>  -  <y,Ly>  (15) 

The  above  expression  can  be  integrated  in  two  different  ways  and  can  also 
be  written  in  terras  of  boundary  conditions  and  initial  conditions.  It  is 
assumed  that  these  boundary  conditions  are  assigned  in  such  a  manner  that  th 


^Stacey,  W.  M.,  Jr,  Variational  Methods  in  Nuclear  Reactor  Physics,  Academic 
Press,  1974. 


4 


above  bilinear  concomitant  is  identically  zero  for  all  independent  variables, 


i  .e . , 


D  =  0 


(16) 


Then  the  first  variations  of  D  also  vanish. 

*  SD  =  SD(6y)  +  6D(6y)  =  0  (17) 

Since  <5y  and  6y  are  independent  of  each  other,  then 

SD(6y)  =  <Sy,Ly>  -  <y,L6y>  =  0  (18) 

SD(6y)  =  < y , L 6 y>  -  <6y,Ly>  =  0  (19) 

Equation  (19)  is  identical  to  Eq.  (9),  which  is  the  assumed  equality 
previously.  The  implication  is  that  if  Eq.  (16)  is  true  then  Eq.  (9)  or  (19) 
is  automatically  true. 

Since  Eq.  (15)  can  be  expressed  in  terms  of  some  integrals  involving 
boundary  conditions,  Eq.  (16)  can  be  true  if  these  boundary  conditions  are 
satisfied.  The  next  section  will  discuss  integral  of  bilinear  expression  and 
its  boundary  conditions. 


INTEGRAL  OF  BILINEAR  EXPRESSION 

The  integral  of  a  bilinear  expression  for  a  two-dimensional  second  order 
problem  in  space-time  can  be  written  as 

xb  tb 

I  =  /  /  #y(x,t)  ,y(x,t)  ]dt  dx  (20) 

xo  to 

where  [ y , y ]  is  a  given  bilinear  expression  in  the  form 

'Ky.y]  =  aytyt  +  Byty  +  yyyt  +  ^YxYx  +  v>yxy  +  vyyx  +  £yy  (21) 
The  subscripts  t  and  x  indicate  the  partial  derivatives  of  the  function  y  and 

y* 
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Equation  (20)  can  be  integrated  by  parts.  Two  different  forms  of  inte¬ 
gration  and  end  conditions  can  be  obtained.  The  first  form  of  the  integral  is 

xb  tb  _  xb  -  tb  tb  -  xb 

I  =  -J  /  yLydtdx  +  /  (ayt+yy)y|  dx  +  J  Uyx+^y)y|  dt  (22> 

XQ  tQ  Xq  tQ  tQ  Xq 

which  is  obtained  by  integrating  by  parts  on  the  adjoint  variable.  On  the 
other  hand,  we  can  perform  integration  on  the  original  variables  to  give 

xb  ^b  —  xb  —  —  tb  . tb  —  —  .  xb 

I  =  -/  /  yLydtdx  +  /  (ayt+3y)y|  dx  +  J  Uyx+uy)y|  dt  (23) 

Xq  tQ  Xq  tQ  t0  Xq 

where 

Ly  -  (ayt)t  -  Syt  +  (Yy)t  +  (£Yx)x  "  ^Yx  +  (vy)x  “  (2^) 

and 

Ly  =  (ayt)t  +  (3y)t  -  YYt  +  (^yx)x  +  (^y)x  “  vYx  “  eY  (25) 

For  a  two-dimensional  second  order  system  in  space— time  domain,  Eq.  (13) 
becomes 

xb  t-b  -  xb  *-b  — 

D  =  /  /  yLydtdx  -  J  J  yLydtdx  (26) 

Xq  tQ  Xq  tQ 

By  equating  Eqs.  (22)  and  (23)  and  solving  for  D  in  Eq.  (26),  we  are  convert¬ 
ing  the  double  integral  into  two  simple  integrals  in  terms  of  the  boundary 
conditions . 

We  can  express  the  quantity  D  as  the  sum  of  two  parts  Dj  and  D2  as 

D  =  Di  +  D2 
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The  terms  in  involve  the  initial  conditions  of  y  and  y  as 

xb  -  -  -  - 

Di  =  /  {«b(ytbyb-ytbYb)  ~  «o(ytoyo-ytoyo) 

xo 

+  (Tb_i3b)ybyb  ~  (y0-30)yoyo}dx  (27)- 

The  terms  in  D2  involve  the  boundary  conditions  of  y  and  y  as 

tb  -  -  -  - 

D2  =  /  Ub(yxbyb-yxbYb)  -  ^0(yxoyo-yxoyo) 
to 

+  (vb-tib)ybyb  -  ( vo-^o)yoyo^dt  (2g) 

In  order  that  D  =  0  in  Eq.  (16),  it  requires  that 

Di  =  0  (29a) 

and 

D2  =  0  (29b) 


END  CONDITIONS  FOR  THE  ADJOINT  SYSTEMS 

We  may  take  four  different  cases  in  discussing  the  end  conditions  for  the 
adjoint  systems  in  order  to  satisfy  the  requirements  in  Eqs.  (29a)  and  (29b) 
(a).  The  Wave  Equation:  In  this  case  Eq.  (24)  becomes 

Ly  =  (ayt)t  +  (£yx)x  =0  (30) 

and  the  coefficients  are 

Yb  =  ^b  >  Yo  =  ^o>  vb  =  l*b>  vo  =  ho 

cxb  *  0,  aQ  *  0,  £b  ^  0,  and  £0  *  0  (31) 
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Let  us  assume  that  the  adjoint  variables  are 

yb  =  ky0»  y0  =  kyb  (32) 

ytb  =  _ab"laokyto.  yto  =  ~ao~ labkytb  (33) 

yxb  =  ~^b  l£okyxo  an<i  yXo  =  -^o_1^bkyxb  (34) 

where  k  is  constant. 

The  above  boundary  values  satisfy  the  requirement  that  =  D2  =  0  in 
Eqs.  (27)  and  (28).  Thus  it  also  satisfies  Eq.  (16)  that 

D  =  0 

(b) .  Heat  Equation:  In  this  case  Eq.  (24)  is 

Ly  =  -0yt  +  ( yy)t  +  Uyx)x  =  0  (35) 

and  the  coefficients  are 

vb  "  Ub>  vo  =  Vo,  ab  =  0,  Oq  =  0 

Tb  *  3b,  Y0  *  30,  £b  *  0,  £0  *  0  (36) 

Let  the  adjoint  variables  be 

yb  =  (Y0-30)ky0»  y0  =  (Yb~3b)kyb  (37) 

yxb  =  -Ab  l£o(Yb-3b)kyxo»  yxo  =  ~ ^o~ 1 ^b ( To- 30 ) kyxb  (38) 

We  also  have  =  D2  =  0  and  D  =  0. 

(c) .  First  order  partials  of  x  in  Eq.  (24)  are  missing,  i.e., 

vb  =  Vb>  vo  =  Vo>  Yb  *  3b,  Y0  *  30  (39) 

Let 

yb  =  (Y0-30)ky0,  y0  =  ( Yb“3b)ky0  (40) 

ytb  =  ~ab-lcto(yb-3b>kyto>  yto  =  -ao-1%(  Y0-20)kytb  (41) 

yXb  =  -*b-1*o(Tb-3b)kyxo>  yxo  =  -V1  Vy0-30)kyxb  (42) 
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We  also  have  =  D2  =  0  and  D  =  0. 

(d).  First  order  partials  of  t  in  Eq.  (24)  are  missing. 


Yb  =  Bb»  Y0  =  30, 

^b  *  ^b»  vo  *  ^o 

ab  ^  0  >  ao  ^  0  > 

&b  *  0,  and  £0  *  0 

(43) 

Let 

o 

/—S 

o 

1 

o 

Vw' 

II 

1 

y0  =  (vb-Mb>kyb 

(44) 

Ytb  ”  “ab  1(Xo(vb  ^>)kYto> 

yto  =  -ao-1«b(  vo-ho)kytb 

(45) 

Yxb  “  “V1  Vvb“Vib)kYxo> 

yxo  -  -^0-1^b(y0-uo)kyxb 

(46) 

We  also  have 

Di  =  D2  =  0,  and  D  =  0. 

From  the  expression  D  =  0  in  Eq. 

(15)  we  can  conclude  that  6D  =  0 

in  Eq 

(19) 

and  Eq. 

(9)  hold,  which  leads  to 

the  condition  in  Eq.  (10)  that 

6  J 

-  0 

for  all  arbitrary  variations  Sy  and  6y. 

FIRST  VARIATION 

Since  the  variations  Sy  and  Sy  are  independent  of  each  other,  the  part  of 
Sj  in  Eq.  (10)  with  variation  Sy  can  be  expressed  as 

xb  tb  -  xb  tb  - 

Sj(Sy)  =  j  j  SyLydtdx  +  /  J  SyQdtdx  =  0  (47) 

Xq  tQ  Xq  t0 

where  Ly  is  given  in  Eq.  (24)  and  contains  second  order  partial  differentials 
in  y.  It  is  intended  to  include  only  first  order  partial  differentials  and 
the  function  y  itself  in  SJ(Sy).  This  can  be  achieved  by  considering  the 
variation  of  the  bilinear  expression  I  in  Eqs.  (20)  and  (21)  which  gives 

SI  =  Sl(Sy)  +  Sl(Sy)  (48) 
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where 


xb  ,tb  - 

6l(6y)  =  /  /  [(ayt+Yy)6yt  +  (  gyt+eyf  uyx)  <$y  +  ( ^yx+vy) 6yx] dtdx  (49) 


xo  co 


A  different  version  of  the  above  variation  can  be  obtained  from  Eq.  (22)  as 


xb  -  xb  tb 

6I(6y)  =  -J  6yLydtdx  4-  J  5y(  ayt+Yy)  |  dx 
xo  xo  to 

,tb  -  ,  xb 

+  J  6y(  £yx+vy)  |  dt 


(50) 


to  *o 

Equating  Eqs.  (49)  and  (50),  solving  for  the  term  containing  integral  for  6yLy 

and  substituting  into  Eq.  (47)  we  have 

xb  -  tb  tb  -  xb 

6  J(  6  y)  =  /  (ayt+yy)6y|  dx  +  /  Uyx+vy)<$y|  dt 
*o  tQ  tQ  xQ 


SyQdtdx 


xb  tb  - 

-  J  J  [ocyt+Yy)6yt  +  (eyt+ey+*uyx)6y  +  (*yx+vy)fiyxldt:dx  =  0  (51) 

xo  tQ 

This  is  the  key  equation  which  uses  variational  principle  in  solving  a 
mixed  initial  and  boundary  value  problem.  The  above  equation  contains  only 
6y,  6yt>  and  6yx  and  none  of  the  variations  of  higher  derivatives.  The 
dependent  variable  contains  only  y,  yt,  and  yx  and  no  higher  partials. 


TRANSFORMATION  OF  COORDINATES 

The  integral  signs  in  Eq.  (51)  can  be  converted  into  summation  signs  if 
discrete  intervals  for  integration  are  used.  We  may  take  some  scale  factors 
to  nondimensionalize  the  problem  by  giving 
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to  =  0,  tb  =  1  0  <  t  <  1  (52) 

xo  =  0,  xb  =  1  0  <  x  <  1  (53) 

Moreover,  Eq.  (51)  can  be  discretized  by  letting 

C  =  Ht  -  i+1  0  <  £  <  1  i  »  1,2 . H  (54) 

TV  =  Kx  -  j+1  0  <  n  <  1  j  =  1,2 . K  (55) 

where  H  and  K  are  number  of  intervals  for  t  and  x  respectively. 

Thus  the  partial  derivatives  are 

3y  3y 

yt '  Tt ' H  ii  =  Hy«  <56> 

8y  8y 

y*  '  ii  ‘  K  in  '  Ky"  (57) 

Use  of  Eqs .  (52)  through  Eq.  (57)  then  leads  to 


K  1  tb  1 

oJ(6y)  =  )_  /  [c*Hyr(t,j)  +  yy( i , j ) ] $y( i , j )  I  dp 

j=l  0  tQ  K 

H  |  _  ^-j_^  I 

+  l  /  UKy^i.j)  +  vy(i,j)]6y(i.j)  I  -d£ 

i=l  0  Xq  H 

K  i  H  i  i  , 

+  l  /  {  l  /  6yC 1 , j ) q  -  d£}  -  dn 

j  =  l  0  i=1  0  H  K 


K  1  H  1 

~  l  !  (  l  /  [ ( aHy r ( i » J )  +  yy(i» j ) )HSyr(i , j ) 

j  =  l  0  i=1  0  ^ 


+  ( , j )  +  ey( i , j )  +  MKyn( 1 » j  )  )  SyU » j  ) 

+  UKy^i.j)  +  vy(i»j))K6yn(i»j>]  -  dO  -  dn  (58) 

H  H 
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GRID  SYSTEMS 


The  (16x1)  vector  has  a  grid  of  four  (4x1)  vectors  through 

Y^( > j ) ,  thus 

YCi.j)  =  t [Y1(i»J>]T  [Y2(i’j)]T  [Y3C ) ]T  [Y4(i,j)]T  (59) 

Each  of  the  (4x1)  vectors  has  four  components,  consisting  of  the  function,  its 
first  partials  in  both  directions,  and  its  mixed  partial. 


*1 


(i,j)  = 


y(Ci,nj) 

yUi,nj+i) 

y^(Ci,nj) 

Y3(i>j)  = 

y£(Si>hj+i) 

yri^i*11^^ 

yn^i’^j+O 

yCn^i’1^ 

y5n^i»nj+l) 

(i,j)  = 


y(5i+i,nj) 

y(Ci+i,nj+i) 

ye<5i+i,nj) 

Y4(i>j)  = 

y^(5i+i ,nj+i) 

yn^i+i»nj) 

yn^i+l»nj+0 

y^ri^i+l  ,rij^ 

ycn^i+i*nj+i^ 

Y2 


(60) 

If  we  increase  the  row  index  from  i  to  i+1 ,  then  the  grid  point  shifts  down  by 
one  step  and  the  following  holds 

Yl(i+l,j)  =  y2^ 1 > d )  Y3(i+1>j)  =  Y4^ 1 > j  ^  (61) 

If  we  increase  the  column  index  from  j  to  j+1  then  the  grid  point  shifts  to 
the  right  by  one  step  and  one  obtains 

Yl(i,j+1)  =  y3(i,j)  Yz(i,j+1)  =  y4^ 1 > j )  (62) 
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The  following  diagram  shows  the  relationship  of  the  grid  system. 


Yi^.j) 

=  Yi^.j+l) 

Y3^>3+1) 

[y(i.i+l) ] 

Yo( 1 > 3 ) 

YaO,  j  +  D 

Y^O+l  ,3) 

Y3(i+l,j)  =  Y1(i+1>d+1) 

Y3(i+1 ,3  +  1) 

[Y(i+1»3)] 

[y(i+l ,3+1)] 

Y2^i+1  >j) 

Y4^i+l>j)  =  Y2^i+^>j  +  ^ 

Y4^i+1  >3  +  0 

SPLINE  FUNCTION 

We  may  express  the  variables  y(^>j)  and  5y(i>j)  in  Eq.  (58)  in  terms  of 
the  (1x16)  spline  function  aT(£,u)  and  the  (16x1)  node  point  function 
as  follows. 

y(i>j)U,n)  =  aT(5,n)Y(i,j)  (63) 

where 

aT(5,n)  =  {[a1(£,n))T  [a2(5,n)]T  [a3(£,n)]T  [a4(5,n)]T  (64) 

and 

Sy(i»3)(5,n)  =  aT(£,n)SYU,j)  (65) 

A  typical  terra  for  a  product  can  be  written  as 

SyOj  )y(i,j  )  =  [6Y(i>j)]Ta(5,n)aT(5,'l)Y(i»j)  (66) 
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BICUBIC  HERMITE  POLYNOMIAL  SPLINES 


With  the  aid  of  Eq.  (59),  Eq.  (63)  may  be  expressed  as 

y(i»j)(5,n)  =  [a1(5tn)]T^i^1»^)  +  [a2(5,n)]TY2(i’j) 

+  [  a3  (5  >i"l)  ]  tY3^  1 » j  ^  +  [a4(5,n)]TY4^1»j^  (67) 

The  bicubic  Hermite  polynomial  spline  is  continuous  in  the  functional  value, 
its  first  partials  in  two  directions,  and  its  mixed  first  partial  in  both 
directions.  The  bicubic  Hermite  polynomial  spline  gives 


—  — 

T 

<t> ( 5 )  <Kn) 

<K 5)  p(n) 

t[5)  <Kn) 

a3(5,n)  = 

<K5)  p(n) 

<J>(5 )  t(n) 

<t>(5)  “(h) 

^(5)  t(n) 

^(5)  m(n) 

where 


p(5)  <Kn) 

T 

p(5)  p(n) 

u(5)  <t>(n) 

m(5)  p(n) 

p(5)  'J'(n) 

a4(5,n)  = 

p(5)  w(n) 

u(5)  Kn) 

m(5)  m(n) 

_  _ , 

—  — 

(68) 


<f>(5)  =  1  -  352  +  253 
¥(5)  =  5  ~  252  +  53 
p(5)  =  352  ~  253 
w(5)  =  -52  +  53 


♦5(5)  =  "65  +  652 

+5(5)  =  1  -  45  +  352 

P5(t)  =  65  -  652 

m^(5)  =  -25  +  352  (69) 
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At  grid  points  (nodes)  the  value  of  (  or  n  takes  the  value  of  0  or 
have 

(f)(0)  =  1  1|>(0)  =  0  (f)(1)  =  0  1|»(1)  =  0 


<f>£(0)  =  0  ip£(0)  =  1 

p(0)  =  0  u(0)  =  0 

P?(0)  =  0  U5(0)  =  0 


4>c(l)  =  0 

p(l)  =  1 
P£(l)  =  0 


*5<1)  =  o 

u(l)  =  0 
u^l)  =  1 


Thus  we 


(70) 


It  is  noted  that  the  diagonal  elements  of  the  matrix  are  unity  and  the  off 
diagonal  terras  are  zeroes.  Similar  expressions  are  held  for  <{>(ri),  etc.  in 
terras  of  n .  For  example 

<f>(n)  =  l  -  3n2  +  2n3  ,  etc.  (71) 


CONSISTENCY  AT  NODES 

To  show  that  Eqs.  (67)  through  (69)  are  consistent  at  the  node  points,  we 
will  check  only  the  following  cases. 

(1)  For  the  case  5  =  0  and  r|  =  0,  from  Eqs.  (68)  and  (70)  we  have 
a2 (0 ,0)  =  a3(0,0)  =  a4(0,0)  =[0000] 

a1 (0 ,0)  =[1000]  (72) 


(a) 


5=0  =  [al(5,n)]T5=o 
n=0  n=0 


y(i»j)(5,q) 

=  [100  OlYi^.j)  =  y(i,j)(o,0) 


(73) 
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(b) 


(c) 


5=0  =  [a.1^C^,n)  ]  T^=0 
n=0  n=0 


+5(5)  <Kn) 

T 

y(i»j)  (0,0) 

i|>5(5)  <Kn) 

75(1, j ) (0,0) 

<|)5(5)  iKn) 

Yr/MHO.O) 

'['5(5)  iKn) 

5=0 

y5n^i#^(0,0) 

n=0 


=  [010  0]  1 » J  ^ 

=  y^iJHOjO) 


y5n(l,J)<5,n) 


5=0  =  [a15n(5,n)]T5=o  Y^1*^ 
n=0  n=0 


<l>5(5)  4>n ( n ) 

T 

y(i » j ) (0 ,0) 

^5(5)  <|>n(n) 

y5^i>J^(0,0) 

<t>5(5)  '|)n(Ti) 

yn(i,j)(0,0) 

^5(5)  4’n(r') 

5=0 

y5n(i,j)(0,0) 

n=0 


(74) 


=[000  1 ] Y}( 1 > J ) 

=  y5n^i,:’)(0>0) 

The  above  expressions  show  that  the  function,  its  first  partial  in  one 


(75) 


direction,  and  its  mixed  partial  are  consistent  and  continuous  at  the  node 
point  5=0  and  n  =  0. 
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(2)  For  the  case  5  =  1  and  n  =  1,  from  Eqs.  (68)  and  (70)  we  have 
a^d.l)  =  a2(l,l)  =  a3(l,l)  =  [0  0  0  0] 

a4(l,l)  =[1000]  (76) 

(a)  Thus 

y(i»j)(5,n)  £;=!  =  [a4(5,n)]T5=i  Y^.J) 

n=l  n=l 

=  [1  0  0  0]Y4d,j)  =  yd,j)(l,l)  (77) 


(b) 


y^(i»j)(5,n) 


5=1  =  [a45(5,n)]T5=i  Y4(1»j> 
n=l  n=l 


(c) 


p^(5)  p(n) 

T 

y(i*j)(l,l) 

w^(5)  p(n) 

y^Ci.dd.l) 

P^(5)  w(n) 

yn(i>j)(i,i) 

<^(5)  w(5) 

5=1 

Y5nd»j ) ( 1 » 1) 

n-1 


=  [010  0]  Y4^  1 » j  )  =  y^d ,  j  )  (  1 , 1) 


y£n(i,j)(dn) 


5=1  =  [a45nd,n)]T5=i  Y4d,j) 


n=l 

n=l 

p^(5)  Pn(n) 

T 

yd,j)(l,l) 

•^(5)  pn(n) 

y5(i»J)(l,l) 

p^(5)  un(n) 

yn(i.j)(i,i) 

0)^(5)  wn(n) 

5=1 

y5nd.j)(l,l) 

n=l 


=  [0  0  0  l]Y4d,j)  =  y^n(i»j)(l,l) 


(78) 


(79) 
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It  can  be  proved  that  all  the  16  elements  at  four  corners  of  the  grid  are 


consistent.  It  can  also  be  proved  that  the  function,  its  two  directional 
first  derivatives  and  its  mixed  partial  are  continuous  at  all  grid  points. 


WAVE  EQUATION 

Let  us  take  a  special  care  for  further  study  of  the  mixed  initial  and 
boundary  value  problems.  We  choose  the  wave  equation  where  the  parameters  in 
Eq.  (24)  are 

3=y=U=v=e=0  (80) 

and  a  =  const  *  0  (81) 

%  -  const  *  0  (82) 

Then  Eq.  (24)  becomes 

Ly  =  aytt  +  *yxx  =  ~Q  (83) 

Eq.  (58)  is  simplified  to 

-  K  i  -  #  t\y 

6j(6y)  =2,  —  /  6y(^>  J  )y^^>  d  ^dn  | 


j-1 


K  0 


H  o  if  i  _  xy* 

+  l  --  /  6yCl.J)ynCl,j)d5| 


i-1 


H  0 


K  H  i  li- 
+  I  I  —  ff  6y( 1  >3  )qd5dri 
j-1  i-1  HK  0  0 

K  H  ii  _ 

+  I  I  /  /  [ —  6y^n,j)y^(i,j) 


j-1  i-1 


0  0  K 


£K  — 

+  —  Sy^O*  j  )yT1(i»j)  ]  d^dri  =  0 
H 


(84) 
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Differentiating  Eqs«  (63)  and  (65),  and  substituting  into  Eq.  (84)  we 


have 

6j(6y)  =  l  --  1 6Y^ 1 » J )  1 T/L  a(5,n)a?(^,n)dn|  yO.J) 
2=\  K  0  tQ 

H  ^  ^  _  T  1  ^b 

+  y  —  [  SYCi.  j )  ]  /  aU,n)anU,n)d5|  YU,j) 

i=l  HO  xQ 


K  H  i  -  Til 

+  T  T  —  /  /  aU,n)QU,n)d!-dn 

j=l  i=l  0  0 

K  H  rvu  —  Til  /  »  .  \ 

+  l  l  —  [SYO.j)]  /  /  a5(5,n)a?T(5,n)d5dn  Y(i>3) 

n  =  i  -r  =  i  H  0  0 


K  H  pv  _  Til  /  *  •  \ 

+  I  l  —  [6Y(i»J)]J  /  an(5,n)anT(5,n)dWn  Y<i»J>  =  0  (85) 

j=l  i=l  H  00 

Eq.  (85)  may  be  written  into  a  different  form  as 

K  t 

6j(6y)  =  l  [ 6Y( tb » j ) ]  P05(tb)Y^tb»J) 

j-1 

-  I  [6Y(t0»JhTP05(t0)^to»J) 

j=l 

H  x 

+  l  [6Y(i*xb>]  P0n(xb)Y(i»xb) 
i=l 

H  x  .  . 

-  y  [SyOjXq)]  p0r)(Xo)Y(i>xo) 

i=l 

K  H  x 

+  l  l  [6Y<i»J)]  q(i»d) 
j=l  i=l 

K  H  -  T 

-ll  [  6yU > j)  ]  p( i » j ) yO » j  )  =  0  (86) 
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It  is  noted  that  the  first  two  terras  involve  initial  values,  the  next  two 
terras  involve  boundary  values,  and  the  last  two  terras  involve  interior 
quantities  within  the  region. 

Equation  (86)  uses  the  following  notations 

ctH  1 

Poe(tb)  =  —  /  aU,n)a?TU,n)dn|  (87) 

K  0  s  t=tb 


aH  ,1  m  i 

PoeUq)  =  —  J  a(5,n)ari(5,h)dn| 

^  K  0  t=tr 


JlK  1  T 

P0n(xb>  =  —  I  aU,n)an  U,n)d5| 

1  K  0  x=xb 


1 

P0n(xo)  =  "  L  a(5,n)anT(5,n)d5| 


K  0  X=XQ 


aH  £K 

p  =  “  *«  +  "  pnn 


.1  .1 


Pee  =  /  /  arU,n)a5TU,n)d5dn 

J0  0  ^ 


and 


.1  ,1 


Pnn  =  /  /  anU,n)anTU,n)d5dn 


0  0 


q(i,j)  =  --  J1  /  aU,n)QU,n)d£dn 
HK  0  0 


(88) 

(89) 

(90) 

(91) 

(92) 

(93) 

(94) 


For  a  given  spline  function,  such  as  bicubic  Hermite  polynomials,  a(£,n) 
is  given  in  Eqs.  (64),  (68),  and  (69).  For  a  given  grid  the  number  of  node 
points  are  known,  so  are  H  and  K.  Thus  Eqs.  (87)  through  (93)  can  be 
determined  and  stored  in  advance.  For  any  given  forcing  function  Q(5>h),  Eq. 
(94)  can  also  be  evaluated. 
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The  (16x1)  vector  Y^»j )  in  Eq.  (86)  was  defined  in  Eq.  (59).  Its 
components  can  be  overlapped  as  given  in  Eqs.  (61)  and  (62).  We  have  found 
previously  that  the  first  term  in  Eq.  (86)  can  be  dropped  because  it  can 
automatically  satisfy  the  final  conditions  for  an  initial  value  problem.  In 
the  second  term  the  function  is  known  because  the  coefficients  are  the 

initial  values.  Although  some  of  the  boundary  values  Y^*xb^  amd  y( i>xo)  are 
given,  most  of  these  terras  are  to  be  determined.  The  entire  problem  is  to 
solve  for  Y'*- » j  )  by  setting  to  zero  the  assembly  coefficients  of  the 
individual  elements  of  SY^jj). 

It  is  a  tedious  task  to  assemble  these  coefficients.  These  will  be 
performed  in  the  future  as  a  separate  report. 

CONCLUSION 

A  bilinear  form  of  the  original  and  adjoint  variable  is  employed  in 
determining  the  coefficients  of  the  variations  of  the  functions  and  their 
first  derivatives.  There  is  no  terra  involving  the  variations  of  any  higher 
derivatives.  A  bicubic  Herraite  polynomial  is  used  which  gives  continuity  in 
the  functions  and  first  partial  derivatives  in  space  0£  time,  together  with 
the  mixed  first  partial  derivative  in  space  and  time.  In  solving  mixed 
boundary  and  initial  value  problems  of  a  second  order  partial  differential 
equation  using  spline  functions,  the  computation  may  be  simplified 
considerably  if  the  variable  in  time  can  be  truncated  into  arbitrary  sections. 
The  entire  problem  is  divided  into  several  strips  of  distinct  time  intervals, 
each  strip  containing  mostly  the  boundary  value  problem. 
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The  variational  principle  for  spatial  and  temporal  problems  with  boundary 
and  initial  conditions  has  been  investigated.  This  variational  principle  is 
very  general  in  scope  and  can  be  applied  to  many  linear  partial  differential 
equations.  The  principle  is  applicable  if  the  bilinear  concomitant  is 
identically  zero.  This  leads  to  the  requirement  that  a  set  of  end  conditions 
for  the  adjoint  systems  must  be  found  to  satisfy  this  condition.  Otherwise 
the  variational  principle  as  stated  may  not  be  applicable. 

Both  the  wave  equation  and  the  heat  equation  (with  one  dimensional 
spatial  direction)  satisfy  these  variational  principles.  For  future  work  the 
analytic  solution  of  these  equations  using  the  finite  element  method  will  be 
studied.  The  assembly  of  the  elements  of  the  matrices  involved  in  the 
formulation  will  be  demonstrated.  The  stability  problem  in  numerical 
solutions  on  these  equations  will  also  be  investigated.  This  lays  the 
foundation  for  the  gun  dynamics  problem  to  be  studied  in  the  future. 
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ATTN:  DOCUMENTS  SECT,  BLDG  4484  2 

REDSTONE  ARSENAL,  AL  35898 


COMMANDER 
REDSTONE  ARSENAL 

ATTN:  DRSMI-RRS  1 

-RSM  1 

ALABAMA  35809 

COMMANDER 

ROCK  ISLAND  ARSENAL 

ATTN:  SARRI-ENM  (MAT  SCI  D1V)  1 

ROCK  ISLAND,  IL  61202 

COMMANDER 

HQ,  US  ARMY  AVN  SCH 

ATTN:  OFC  OF  THE  LIBRARIAN  1 

FT  RUCKER,  ALABAMA  36362 

COMMANDER 

US  ARMY  FGN  SCIENCE  &  TECH  CEN 
ATTN:  DRXST-SD  1 

220  7TH  STREET,  N.F.. 

CHARLOTTESVILLE,  VA  22901 


COMMANDER 

US  ARMY  MATERIALS  &  MECHANICS 
RESEARCH  CENTER 
ATTN:  TECH  LIB  -  DRXMR-PL 

WATERTOWN,  MASS  02172 


1 

(PLASTICS  TECH  1 

EVAL  CEN') 

1 

1 

1 

1 

(STINFO)  2 


NOTE:  PLEASE  NOTIFY  COMMANDER,  ARRADCOM,  ATTN:  BENET  WEAPONS  LABORATORY, 

DRDAR-LCB-TL ,  WATERVLIET  ARSENAL,  WATERVLIET,  N.Y.  12189,  OF  ANY 
REQUIRED  CHANGES. 


TECHNICAL  REPORT  EXTERNAL  DISTRIBUTION  LIST  (CONT.) 


NO.  OF 


COMMANDER 

IJS  ARMY  RESEARCH  OFFICE- 
P.O.  BOX  12;  11 

RESEARCH  TRIANGLE  PARK,  NC  27709 
COMMANDER 

IJS  ARMY  HARRY  DIAMOND  LAB 
ATTN:  TECH  LIB 
2 BOO  POWDER  MILL  ROAD 
ADELPHIA,  ME  20783 

DIRECTOR 

US  ARMY  INDUSTRIAL  BASE  ENG  ACT 

ATTN:  DRXPI-MT 

RCJCK  ISLAND,  IL  61201 

CHIEF,  MATERIALS  BRANCH 
US  ARLwOf  R&5  GROUP,  EUR 
BCK  65,  FPO  N.Y.  09510 

COMMANDER 

NAVAL  SURFACE  WEAPONS  CEN 
ATTN:  CHIEF,  MAT  SCIENCE  DIV 
DAHLGREN,  VA  22448 

DIRECTOR 

US  NAVAL  RESEARCH  LAB 
ATTN:  DIR,  LECH  DIV 

CCTE  26-27  (DOC  LIB) 
WASHINGTON,  3.  C.  20375 

NASA  SCIENTIFIC  &  TECH  INFO  FAC. 

P.  0.  BCK  3757,  ATTN:  ACQ.  BR 
BALTIMORE/WASHINGTCN  INTL  AIRPORT 
MARYLAND  21 240 


COMMANDER 

DEFENSE  TECHNICAL  INFO  CENTER 
1  ATTN:  DTIA-TCA 

CAMERON  STATION 
ALEXANDRIA,  VA  22314 

METALS  &  CERAMICS  INFO  CEN 

1  BATTELLE  COLUMBUS  LAB 

505  KING  AVE 
COLUMBUS,  OHIO  43301 

MECHANICAL  PROPERTIES  DATA  CTR 
BATTELLE  COLUMBUS  LAB 

1  505  KING  AVE 

COLUMBUS,  OHIO  43201 

MATERIEL  SYSTEMS  ANALYSIS  ACTV 
1  ATTN:  DRXSY-MP 

ABERDEEN  PROVING  GROUND 
MARYLAND  21005 


1 


1 

1 


1 


NOTE:  PLEASE  NOTIFY  COMMANDER,  ARRADCCM,  ATTN:  BENET  WEAPONS  LABORATORY, 
DRTAF.-LCB-TL,  WATERVLIET  ARSENAL,  WATERVLIET,  N.Y.  12189,  OF  ANY 
REQUIRED  CHANGES. 


